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ABSTRACT 

We study the dynamical evolution of globular clusters containing populations of primordial binaries, 
using our newly updated Monte Carlo cluster evolution code with the inclusion of direct integration 
of binary scattering interactions. We describe the modifications we have made to the code, as well as 
improvements we have made to the core Monte Carlo method. We present several test calculations 
to verify the validity of the new code, and perform many comparisons with previous analytical and 
numerical work in the literature. We simulate the evolution of a large grid of models, with a wide 
range of initial cluster profiles, and with binary fractions ranging from to 1, and compare with 
observations of Galactic globular clusters. We find that our code yields very good agreement with 
direct A'^-body simulations of clusters with primordial binaries, but yields some results that differ 
significantly from other approximate methods. Notably, the direct integration of binary interactions 
reduces their energy generation rate relative to the simple recipes used in Paper III, and yields smaller 
core radii. Our results for the structural parameters of clusters during the binary-burning phase are 
now in the tail of the range of parameters for observed clusters, implying that either clusters are born 
significantly more or less centrally concentrated than has been previously considered, or that there 
are additional physical processes beyond two-body relaxation and binary interactions that affect the 
structural characteristics of clusters. 

Subject headings: globular clusters: general — methods: numerical — stellar dynamics 



1. INTRODUCTION 

Observations (e.g.. Cote et al. 1996; Bellazzini et al. 
2002b; Rubenstein & Bailyn 1997; Cool & Bolton 2002; 
Bellazzini et al. 2002a), in combination with recent theo- 
retical work (Ivanova et al. 2005), suggest that although 
the currently observed binary fractions in the cores of 
globular clusters may be small (< 10%), the initial clus- 
ter binary fraction may have been significantly larger 
(> 50%). As has been understood theoretically for some 
time, primordial binaries in star clusters generally act 
as an energy source (through super-elastic scattering en- 
counters), producing energy in the core and postpon- 
ing deep core collapse in a quasi-steady state "binary- 
burning" phase. This is analogous to the long-lived main- 
sequence in stars, in which hydrogen is burned to prevent 
collapse. An initial binary fraction of a few percent is 
enough to postpone deep core collapse for many initial 
relaxation times (see Fregeau et al. 2003, for discussion 
and references). In addition to playing a large part in 
the global evolution of a star cluster, dynamical inter- 
actions of binaries also strongly affect the formation and 
evolution of stellar and binary exotica, which include low- 
mass X-ray binaries, recycled pulsars, cataclysmic vari- 
ables, and blue stragglers (Hut et al. 1991; Sigurdsson & 
Phinney 1995; Davies & Hansen 1998; Rasio et al. 2000; 
Ivanova et al. 2005). 

Similarly to dynamical binary interactions, physical 
stellar collisions (both direct star-star collisions, and 
those mediated by resonant binary interactions) also play 
an important role in the evolution of globular cluster pop- 
ulations. Stellar collisions are thought to be one of the 
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two primary mechanisms by which blue stragglers are 
created in dense star clusters (e.g., Mapelli et al. 2004). 
A runaway sequence of stellar collisions of massive main 
sequence stars early in the lifetime of a dense cluster 
may yield a very massive star (> 10'^ Mq), which may 
then become an intermediate mass black hole (see, e.g., 
Giirkan et al. 2006, for discussion and references). 

The complete picture of star cluster evolution can be 
rather complicated, as it includes, in addition to the two 
physical processes just mentioned, single star evolution, 
binary star evolution, and tidal stripping due to the field 
of the Galaxy. Even if one ignores these additional pro- 
cesses, modeling a dense stellar cluster with primordial 
binaries still presents a formidable computational chal- 
lenge. There are at least two reasons for this: 1) dynami- 
cal interactions of binaries — which are typically resonant, 
lasting for many orbits — must be resolved on their natu- 
ral timescale, which is orders of magnitude shorter than 
the cluster relaxation time, and 2) primordial binaries 
extend the life of a cluster. Although the GRAPE se- 
ries of special-purpose computers is steadily increasing 
in performance, direct A'^-body simulation of the evolu- 
tion of clusters with more than a few percent binaries and 
a moderate number of stars (^ 10^) is still quite compu- 
tationally expensive, with computational timescales on 
the order of months. Faster, more approximate meth- 
ods, such as the anisotropic gas model or direct solution 
of the Fokker-Planck equation, suffer from the difficulty 
inherent in incorporating physics beyond two-body re- 
laxation, such as stellar evolution or binary interactions, 
in these methods. The Monte Carlo method bridges the 
gap between these two computational extremes, since it 
allows for the relatively facile inclusion of additional lay- 
ers of physics, provides a star-by-star description of a 
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cluster, and is computationally inexpensive. 

In previous studies using approximate methods like 
Monte Carlo or Fokker-Planck, binary interactions had 
generally been treated using recipes culled from the re- 
sults of large numbers of numerical scattering experi- 
ments. (The work of Giersz & Spurzem (2003), which 
incorporates direct integration of binary interactions, is 
one notable exception.) The recipes are typically known 
only for equal-mass binary interactions, thus prohibit- 
ing the use of a cluster mass function. Thus in order to 
model realistic clusters, which contain a wide range of 
masses, one must numerically integrate each binary in- 
teraction in order to resolve it properly. We have now 
incorporated into our Monte Carlo code a dynamical in- 
tegrator to exactly integrate dynamical interactions of 
binaries, allowing us to evolve clusters with mass spec- 
tra, and perform more realistic comparisons with direct 
iV-body calculations. As we demonstrate below, the di- 
rect integration of binary interactions reduces their en- 
ergy generation rate relative to the simple recipes used 
in Paper III, and yields smaller core radii. 

This is the fourth paper in a series studying the evolu- 
tion of globular clusters using the Monte Carlo method. 
Paper I describes the core method and presents several 
test calculations exhibiting the validity of the code ( Joshi 
et al. 2000). Paper II treats the evolution of tidally- 
truncated clusters with mass spectra (Joshi et al. 2001). 
Paper III adds a recipes-based treatment of binary scat- 
tering interactions and considers the evolution of an en- 
semble of clusters of varying initial binary fraction and 
central concentration, finding that even a small fraction 
of binaries in a cluster is sufficient to support the core 
against collapse significantly beyond the normal core- 
collapse time (Fregeau et al. 2003, hereafter Paper III). 
In this, the fourth paper in the series, we describe our 
new code, perform several tests to ensure its validity, and 
perform a large set of simulations of clusters with primor- 
dial binaries, which we compare with previous results in 
the literature and with observations of Galactic globular 
clusters. Section 2 describes our new code in detail, in- 
cluding the additional physical processes we have added 
(numerical integration of binary scattering interactions, 
and star-star physical collisions), as well as the improve- 
ments we have made to the core Monte Carlo method. 
Section 3 presents a few example results, and compares 
them with semi-analytical theory and previous numeri- 
cal calculations. Section 4 describes the trends evident 
in the grid of cluster models we simulated, and compares 
our results with observations. Finally, in section 5 we 
summarize and conclude. 

2. METHOD 

Here we describe in detail our implementation of the 
Monte Carlo numerical method for simulating the evolu- 
tion of dense star clusters. It incorporates many phys- 
ical processes of relevance in dense star clusters, in- 
cluding two-body relaxation, direct physical stellar col- 
lisions, and dynamical interactions of binaries. Each of 
these physical processes is treated sufficiently accurately 
so as to allow for a rather wide mass spectrum (e.g., 
-Wmax/^jfmin ~ 10'^ for a Salpctcr mass function). For 
now we neglect stellar evolution (both single and binary), 
but plan to include it in our code in the near future. 
Those readers uninterested in the technical details of our 



numerical method can safely skip ahead to the next sec- 
tion. 

2.1. Units 

Before describing our method in detail, we discuss a 
necessary formality. We use the standard iV-body system 
of units (Heggie & Mathieu 1986; Heggie & Hut 2003). 
For reasons of convenience, we use two units of time in 
our code. In addition to the standard iV-body unit of 
time (which is roughly the crossing time), we use the 
relaxation time^. Thus our full (over-specified) system 
of units is given by the standard formulae: 
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where Mo is the initial total mass of the cluster, Eo is 
the initial total energy of the cluster, Nq is the initial 
number of stars, and 7iVo is the Coulomb logarithm. The 
quantity 7 is a function of the initial structure of the 
cluster, and is only needed when converting time in code 
units to physical units. Thus it does not need to be 
specified for purely relaxation calculations, which can be 
quoted in units of the relaxation time, but it does need 
to be specified for calculations which include additional 
physics, like physical collisions and binary interactions. 
For our simulations we set 7 via comparisons with iV- 
body results (where available), as discussed later. 

2.2. Standard Definitions 

It is often useful to be explicit about how certain de- 
rived measurable quantities are calculated. For the half- 
mass relaxation time we adopt the standard definition 
(Spitzer & Hart 1971): 
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where N is the number of bound cluster objects (single 
star or binary), is the radius containing half the mass 
of the cluster, and M is the total cluster mass. Most 
plots presented in this paper use as time unit the initial 
half-mass relaxation time, i.e. trh evaluated at the start 
of the simulation. For the core radius we use the density- 
weighted average described in Casertano & Hut (1985), 
with J = 6 and where the averaging is performed from the 
cluster center out to the half mass radius. In this scheme, 
a local density is estimated for each star in the cluster 
by taking the average density within a sphere centered 
on the star with radius at the jth star. The core radius 
is then estimated by taking the local density weighted 
average of star position out to the half mass radius. We 
use the same density-weighted averaging scheme to es- 
timate the central density, velocity dispersion, and av- 
erage mass, except in places where otherwise noted. In 

^ Using the relaxation time as tlic time unit removes A'^ from any 
equations of relaxational evolution. 
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sonic cases, as described in the text, we also plot the core 
radius as measured by the standard definition (Spitzcr 
1987): 
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where ctq is the central three-dimensional velocity dis- 
persion, and po is the central mass density. Note that 
the results we compare with below from Heggie et al. 
(2006) adopt the density-weighted average definition of 
the core radius. As shown in Fig. 13, the two definitions 
differ minimally when the cluster is in the binary burning 
phase. 

2.3. Two-Body Relaxation 

Two-body relaxation is the primary physical process 
responsible for the diffusion of energy in a star cluster, 
and thus for its global evolution (Heggie & Hut 2003). 
We use the Henon orbit-averaged Monte Carlo method 
to simulate two-body relaxation (Henon 1971). For a 
detailed description of the basic method we employ, see 
Joshi et al. (2000). In summary, a timestep in the code 
consists of the following: 

1. Using each star's radial position r and mass m, the 
potential <I>(r) is calculated under the assumption 
of spherical symmetry (each star is represented by 
an infinitesimally thin spherical shell). 

2. Each pair of stars neighboring in radius undergoes a 
hyperbolic "super" -encounter, with scattering an- 
gle chosen so as to represent the cumulative effect 
on each star of many long-range, small-angle two- 
body scattering encounters with all other stars in 
the system. 

3. Using the new radial velocity Vr and tangential ve- 
locity vti the new specific energy E and angular 
momentum J of each star is calculated (using $(r) 
from step 1). 

4. A new position and corresponding velocity is cho- 
sen for each star by picking a point on its orbit 
randomly, sampled in accordance with the amount 
of time spent at each radial position (i.e., weighted 
by 1/vr). 

We have made two improvements to the fundamen- 
tal method which were necessary to accurately treat star 
clusters with wide mass spectra, and for the stability of 
the long-term evolutions needed for clusters with primor- 
dial binaries. 

The first improvement is a rather simple one that pro- 
vides for self-consistency in step 4 above. Solving for 
the new position of a star along its orbit appears to be 
a straightforward matter: using the potential calculated 
in step 1, one writes down the energy equation for the 
orbit, E = <I>(r) + J^/2r^ -|- ^w^, solves for the pericen- 
ter and apocenter, then samples the radial position with 
a weighting inversely proportional to Vr- However, the 
potential calculated in step 1 includes the contribution 
from the star whose orbit we are trying to solve. In other 
words, the star on its orbit feels the gravitational effect 
of itself at its old position. This inconsistency is ignored 
in the standard Monte Carlo method (Henon 1971; Joshi 



et al. 2000), although it is corrected in the new Monte 
Carlo code of Freitag & Benz (2001). Neglecting this in- 
consistency for the case of equal-mass clusters produces a 
minimal effect on the overall evolution, slightly postpon- 
ing core collapse, and leading to a steady drift in the to- 
tal system energy. However, when one considers clusters 
with even modestly wide mass spectra (~ 0.1-10 M©), 
the errors are much less benign. Since the most mas- 
sive stars in the mass spectrum contribute proportion- 
ately more to the cluster potential, it is the calculation 
of their orbits that is most inaccurate. In our simulations 
we have seen that as the mass spectrum is widened, the 
core collapse time gets progressively longer than what 
would be expected — from direct iV-body results and the 
Monte Carlo calculations of Freitag et al. (2006b) — until 
it is prevented completely. In other words, this inconsis- 
tency acts as a spurious energy source, postponing core 
collapse. 

Correcting the inconsistency in the potential is 
straightforward. We simply add a correction term to the 
potential when solving for a star's orbit. For star j with 
mass rrij , originally at position rj when the potential was 
last calculated, the correction term is 



r > Ti 



r < To 
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In principle, for total consistency, one could also add the 
self-gravity of the star, —Gmj/2r, since it is treated as 
a spherical shell. However, we find that adding such a 
term leads to unphysical behavior for clusters with wide 
mass spectra when the orbit of one of the more massive 
stars lies within the innermost few stars (and thus the 
approximation of spherical symmetry breaks down), the 
result being that the star acts as an energy source, ulti- 
mately preventing core collapse. For narrow mass func- 
tions, the addition of the self-gravity term has no notice- 
able effect on the evolution, as found by Freitag & Benz 
(2001). Note that for the results presented in Freitag 
et al. (2006b) and Freitag et al. (2006a), which consider 
core collapse for wide mass spectra, the self-gravity term 
is not included (Freitag 2005, private communication). 

The second improvement to the code concerns energy 
conservation and the long-term stability of the code. 
From the description of a timestep above, it is clear that 
the potential used to find the new positions of stars in 
their orbits (in step 4) lags behind by a timestep. The 
result is a steady drift in the total system energy. This 
can be compensated for by a technique that considers the 
mechanical work done by the potential (since it is chang- 
ing with time) on each star in the system, and uses it to 
more accurately calculate the velocities at the new po- 
sition on the orbit (see Stodolkiewicz 1982, for details). 
Briefly, in this method the new specific kinetic energy 
of each star at its updated position in step 4, ■^v'^^^, is 
augmented by a term that corresponds to the mechanical 
work done by the potential on the star. The tangential 
component of the new velocity is set according to angular 
momentum conservation as wt.now = J/fnew, where J is 
evaluated in step 3. This is the same way in which it is 
set in the standard Monte Carlo method. The radial ve- 
locity is then simply Wr,ncw = (^^now~''^t now)"'^^^- Note the 
asymmetry in which the components of the new velocity 
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are set. Clearly only the radial component of the veloc- 
ity takes into account the mechanical work done by the 
potential. Thus it can happen that v^^^ < new yield- 
ing a nonsensical result for Vi^new, in which case ad hoc 
prescriptions must be used, sometimes leading to spuri- 
ous results. We use a modified version of this method in 
which we simply preserve the ratio Vr.now/wt,ncw as pre- 
dicted by the standard Monte Carlo method and scale the 
velocities so that now + new = '^new This technique 
appears to violate angular momentum conservation for 
individual orbits, but since the Monte Carlo method as- 
sumes spherical symmetry, the total angular momentum 
of the system remains statistically consistent with zero. 
Our modification to the method of Stodolkicwicz (1982) 
yields results for clusters with mass spectra that are more 
consistent with direct A^-body simulations, and the sim- 
ulations of Preitag ct al. (2006b). Moreover, the method 
provides for improved energy conservation throughout 
long cluster runs, typically conserving energy to within 
a part in ^ 10^ over tens of half-mass relaxation times. 

To demonstrate the validity of our improved technique 
for two-body relaxation, we have compared calculations 
of clusters subject only to relaxational evolution with the 
results from other numerical techniques. Fig. 1 shows the 
evolution of the Lagrange radii for a single-component 
Plummer model (model Tl in Table 1) calculated with 
our Monte Carlo code (solid lines), and compared with 
a direct A''-body code (dotted lines). The agreement be- 
tween the Monte Carlo method and direct A^-body is 
clearly excellent for this model. The core collapse time 
of tcc/tih = 17.6 is in good agreement with most other 
approximate techniques, as can be seen in the table of 
Freitag & Benz (2001). For the comparison we had to 
convert the dynamical time units of the A'^-body model 
to relaxation time units, using a value of 7 = 0.10 in 
the Coulomb logarithm. This value is in good agreement 
with theoretical arguments (Henon 1975) and other nu- 
merical calculations (Giersz & Heggie 1994; Freitag & 
Benz 2001). 

Wc have also looked at the evolution of the density 
profile, as shown in Fig. 2, initially and at the time of 
core collapse. A power-law density profile with p cx r~^'^ 
clearly develops at late times. The power-law index of 
—2.3 is in good agreement not only with the results of 
other Monte Carlo calculations (Freitag & Benz 2001), 
but also with A^-body simulations which give an index 
of —2.26 (Baumgardt et al. 2003), and self-similar an- 
alytical Fokkcr-Planck calculations and coarse dynamic 
renormalization calculations which give a power-law in- 
dex of -2.23 (Szell et al. 2005). 

Next we considered the evolution of a model with a 
moderately wide mass spectrum. Fig. 3 compares the 
evolution of the Lagrange radii calculated with our code 
with a direct A^-body calculation of a Plummer model 
with a Kroupa initial mass fimction from 0.1 Mq to 
10 Mq (model T2). Our model used A^ = 10'^ stars while 
the A^-body model used A^ = 131072. Again, the agree- 
ment is quite good, at least to within the level of noise in 
the A^-body simulation. The A^-body model does not un- 
dergo deep collapse, since it appears to form a three-body 
binary that stalls core collapse. (Note that three-body bi- 
nary formation is not included in our Monte Carlo code.) 
A value of 7 = 0.05 in the Coulomb logarithm was used 
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Fig. 1. — Evolution of the Lagrange radii for a single-component 
Plummer model (model Tl in Table 1) calculated with our Monte 
Carlo code (solid lines), and compared with a direct A^-body cal- 
culation (dotted lines). The Lagrange radii shown enclose a fixed 
fraction of the total bound cluster mass of (from bottom to top) 
0.01, 0.02, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, and 0.9. Our 
model had 5 X 10* stars, while the A'^-body model had 65536. The 
time unit in the A'^-body model was converted from dynamical times 
to relaxation times using a value of 7 = 0.10 in the Coulomb loga- 
rithm. 




Fig. 2. — Three-dimensional mass density profiles of model Tl 
initially and at core collapse. At core collapse the cluster develops 
a large power-law profile with exponent —2.3, in good agreement 
with other analytical and numerical calculations. 

to convert between dynamical time units and relaxation 
time units. Fig. 4 shows the evolution of the average mass 
within the Lagrange radii, compared with A^-body. It ap- 
pears that for the innermost Lagrange radii (which have 
the largest average mass) , our Monte Carlo method pre- 
dicts an evolution that lags behind the A^-body method, 
but eventually catches up at late times. Note that the 
Monte Carlo method of Freitag ct al. (2006b) suffers from 
the same malady, albeit to a lesser degree. 

Finally, we considered the evolution of a model with a 
very wide mass spectrum. Fig. 5 compares the Lagrange 
radii with A^-body for a Plummer model with a Salpeter 
initial mass function from 0.2 M© to 120 M© (model T3). 
Our model used A^ = 1.25 x 10^ stars while the A^-body 
model used A = 262144. Again, the agreement is quite 
good, although the A-body model is rather noisy, cs- 
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Fig. 3. — Evolution of the Lagrange radii for a Plummor model 
with a Kroupa initial mass function from 0.1 Mq to 10 Mq (model 
T2) calculated with our Monte Carlo code (solid lines), and com- 
pared with direct A'-body (dotted lines). The Lagrange radii shown 
enclose a fixed fraction of the total bound cluster mass of (from bot- 
tom to top) 0.01, 0.02, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, and 
0.9. Our model had 10^ stars, while the A^-body model had 131072. 
The conversion between dynamical time units and relaxation time 
units required a value of 7 = 0.05 in the Coulomb logarithm. 
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Fig. 5. — The evolution of the Lagrange radii for a Plummer 
model with a Salpeter initial mass function from 0.2 4-/0 to 120 Mq 
(model T3) calculated with our Monte Carlo code (solid lines), 
compared with A'^-body (dotted linos). The Lagrange radii shown 
enclose a fixed fraction of the total bound cluster mass of (from 
bottom to top) 0.01, 0.02, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 
and 0.9. Our model had 1.25 X 10® stars, while the A^-body model 
had 262144. A value of 7 = 0.01 in the Coulomb logarithm was 
used to convert time units. 
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Fig. 4. — Evolution of the average mass within each Lagrange ra- 
dius for model T2 calculated with our code (solid lines), compared 
with Af-body (dotted lines). The values of the Lagrange radii are 
the same as in the previous figure. 
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Fig. 6. — Evolution of the average mass within each Lagrange ra- 
dius for model T3 calculated with our code (solid lines), compared 
with Af-body (dotted lines). The values of the Lagrange radii are 
the same as in the previous figure. 



pecially at late times. Here a value of 7 = 0.01 in the 
Coulomb logarithm was used, which is in good agreement 
with the comparison between Monte Carlo and iV-body 
of Freitag et al. (2006b). Fig. 6 shows the evolution of 
the average mass within the Lagrange radii for the same 
model. Here our model lags even further behind the A^- 
body model at early times, but again catches up at late 
times. Note that the degree to which our model disagrees 
with TV-body appears to be similar to that of the Monte 
Carlo code of Freitag et al. (2006b). 

In general, for two-body relaxation, the agreement of 
our code with the results of direct TV-body calculations, 
as well as those of other approximate techniques, has 
been improved greatly by the two major code modifica- 
tions we have described here. 



2.4. Strong Interactions 

Owing to the flexibility of the Monte Carlo method, it 
is reasonably simple to layer additional physics on top of 
the basic two-body relaxation technique. This includes 
those two-body processes which we call "strong interac- 
tions," including single-single physical collisions, dynam- 
ical binary interactions (between two binaries or a binary 
and a single star), and large-angle scattering. At present 
we have included single-single collisions and binary in- 
teractions. For both we sample the interactions using 
the technique of Freitag & Benz (2002) (also discussed 
in Giersz 2001). In brief, this amounts to evaluating the 
quantity 

^strong — '^♦^^oo'S'itrong'^^ (8) 
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for each pair of stars that neighbor in radial position, 

(12) 

where ^strong the probabiUty for a strong interaction to 
occur, n* is the focal number density of stars (binary or 
single), Woo is the relative velocity of the pair at infinity, 

'5'strong is the cross section for the strong interaction, St 
is the timestep, and the notation "(12)" signifies that 
the first star is of type "1" while the second is of type 
"2" . This quantity is a standard "nav" estimate for the 
interaction probability between stars of type "1" and "2" . 
However, it is the total number density of stars that 
appears in the equation, and not ni or n2. As shown 
in Freitag & Benz (2002), when used in this context, 
eq. (8) yields the correct sampling of the collision rate, 
since the act of choosing the neighboring star samples 
the local density of that type of star (for details, please 
see section 2.4.2 of Freitag & Benz 2002). 
At each timestep, for each pair of stars, the value of 

-^strong i^ evaluated (Pbb for binary-binary interactions, 
Pbs for binary-single, and PcoU for single-single colli- 
sions) and compared with a uniform deviate in X G [0, 1). 
If X < Pstrong the strong interaction is performed, oth- 
erwise the pair undergoes two-body relaxation. 

(12) 

Note that •S'sti-ong can be written in a very general way 
in terms of the maximum value of the classical pericenter 
distance between the pair in their hyperbolic orbit which 
yields a strong interaction. In this case it is 

2GM\ 
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c.(12) 
strong 



2 

max 



1 



(9) 



where fomax is the impact parameter leading to a classical 
pericenter distance of Vp, and M is the total mass of the 
pair. We will use this expression below. 

2.4.1. Single-Single Collisions 

Although the simulations we present in this paper do 
not include physical stellar collisions (wc defer such sim- 
ulations to a future paper), we still include here for 
the sake of completeness a description of the imple- 
mentation of collisions in the code. For direct physi- 
cal collisions between main-sequence stars, the outcome 
can vary greatly depending on Voc- For relative speeds 
greater than the escape speed from the surface of the 
star {voo ^ Vesc ~ 500km/s for a typical solar-mass MS 
star), which can occur in galactic nuclei, a collision typ- 
ically results in a large fraction of the total mass lost 
from the system (see, e.g., Freitag & Benz 2005, and ref- 
erences therein). For v^o ^ Vesc, which is satisfied for 
globular clusters, the result is typically a clean merger, 
with a negligible amount of mass lost from the system 
(e.g., Benz & Hills 1987; Lombard! et al. 2002). Since we 
are concerned with the latter case in this paper, we treat 
physical single-single collisions using the sticky sphere 
approximation, which assumes that if the radii of stars 
touch during a strong interaction they merge with no 
mass loss. The sticky sphere approximation, when used 
for stellar collisions of main-sequence stars in models of 
low velocity dispersion clusters, has been shown to agree 
extremely well with the results of more detailed calcula- 
tions incorporating the results of SPH simulations (Fre- 
itag et al. 2006b, a). With this approximation, the cross 
section for collisions is given by cq. (9) with = R1+R2: 

2GM 
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Fig. 7. — Comparison of the numerically sampled single-single 
star collision rate (circles) with the analytical result (lines) for a 
single-mass Plummcr model with A'^ = 10® stars, for three different 
values of the Safronov number (©o = 1/300, 0.725, and 300). 

where M is the total mass of the two stars. 

We have tested that our code correctly samples single- 
single star collisions so as to reproduce the correct colli- 
sion rate. For a Plummer model, the collision rate can be 
solved for analytically, yielding (Freitag & Benz 2002): 



rf^coll(fl) 

dtdlnR 



-21/4 



X Gn 



i + eo(i + u')i/'l , (11) 



(10) 



where M is the total cluster mass, R is the radial posi- 
tion in the cluster, Rp is the scale radius of the Plum- 
mer model {Rp = 37r/16 in iV-body units), u = R/Rp, 
and &o is the Safronov number (Binney & Trcmaine 
1987). Based on our sticky sphere collision prescrip- 
tion, the Safronov number can be written simply as 
80 = 3Rp/NR^,, where N is the number of stars in the 
cluster, and i?* is a stellar radius. We have extracted the 
collision rate from a simulation of an unevolving Plum- 
mer model (relaxation turned off) composed of = 10^ 
equal-mass stars. Instead of performing collisions, we 
simply recorded them. Fig. 7 shows a comparison of the 
numerically extracted rate (circles) with the analytical 
rate (lines), for three different values of Oq (1/300, 0.725, 
and 300). The agreement is excellent for all values of Qq. 

2.4.2. Binary Interactions 

All dynamical binary interactions (binary-binary and 
binary-single) arc directly numerically integrated with 
Fewbody, an efficient computational toolkit for evolving 
small- A^ dynamical systems (Fregeau et al. 2004). Few- 
body was designed specifically for performing dynamical 
scattering interactions, and thus it is well-suited for our 
purposes. Sec Fregeau et al. (2004) for a detailed descrip- 
tion of the code. Note that Fewbody also properly treats 
physical stellar collisions during binary interactions, us- 
ing the same criterion for a collision as used in the Monte 
Carlo code for single-single collisions. 

For sampling binary interactions, wc use the same tech- 
nique described above for single-single collisions, but 
with the cross sections appropriate to binary interac- 
tions. For binary-binary interactions, the cross section 
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Fig. 8. — Comparison of the numerically sampled binary inter- 
action rates (circles) with the analytical result (lines) for binary- 
binary interactions (black) and binary— single (red) for a single- 
mass Plummer model with N = 10® stars. In this test, binaries 
had the same mass as single stars. This simplified the analytical 
calculation of the rate. 



is given by eq. (9) with rp = Xbb(ao + ai), where Ui 
are the binary semimajor axes, and Xbb is a parame- 
ter. Xbb must be set large enough so that all binary- 
binary interactions of interest are followed. Since in this 
paper we are concerned mainly with the global evolu- 
tion of clusters, we need only follow most of the energy- 
generating binary interactions. In principle, one could 
make Xhh arbitrarily large, so as to capture all poten- 
tially interesting interactions. The result would be many 
more weakly-interacting fly-by interactions, which incur 
an infinitesimal computational cost due to Fewbody's effi- 
cient integration techniques. However, due to the way in 
which the global Monte Carlo timcstcp is chosen (see be- 
low), time would grind to a halt in our code. Clearly, 
then, setting Xbb is a compromise between capturing 
all binary interactions of interest (larger Xbb), and pre- 
venting the timestep from becoming unnaturally small 
(smaller Xbb)- We use Xbb = 2 for the results presented 
in this paper. For binary-single interactions, we take 
Vp ~ XbsO, where a is the binary semimajor axis, and 
set Xbs = 2. We find that the values Xbs = -'fbb ~ 2 
capture almost all the relevant energy-generating binary 
interactions. Test runs with Xbs = ^bb = 4 yield values 
oiTc/rh in the binary burning phase that arc statistically 
consistent with Xbs = ^bb = 2 runs for both small and 
large /b. 

As we did for single-single collisions, wc have per- 
formed a calculation of an unevolving Plummer model 
to test that our code correctly samples the binary in- 
teraction rates. In this test, a fraction fb = 0.5 of the 
N = 10^ stars were binaries. The binaries had the same 
mass as single stars in order to simplify the analytical cal- 
culation of the rate. Again, strong interactions were not 
performed, simply recorded. Since there are two species 
(binaries and single stars), the interaction rate is given 
by eq. (11) with an extra factor of /b for binary-single, 
and for binary-binary. Fig. 8 shows a comparison of 
the numerically sampled binary interaction rates (circles) 
with the analytical result (lines) for binary-binary inter- 
actions (black) and binary-single (red). The agreement 
is excellent. 



Once a binary interaction is deemed to occur, the rel- 
ative velocity at infinity of the pair is taken to be the 
current relative velocity of the pair in the cluster (with 
the angle of each particle's tangential velocity random- 
ized), and the impact parameter, 6, of the interaction is 
chosen uniformly in area out to 6max as given in eq. (9) . 
With all parameters of the scattering encounter set, the 
interaction is numerically integrated with Fewbody until 
an unambiguous outcome is reached. The outcome prod- 
ucts of the interaction are then placed back into the clus- 
ter with their resultant internal properties (mass, binary 
semimajor axis, eccentricity, etc.) and external proper- 
ties (systemic velocity, etc.). 

The only exception to this rule is stable hierarchi- 
cal triples. These stable triples frequently result from 
binary-binary interactions (roughly 20% of the time for 
equal-mass, equal-energy hard binaries; see, e.g., Mikkola 
1983). In principle we could keep the triples in the code, 
and follow their evolution, allowing them to undergo in- 
teractions with other single stars, binaries, and triples 
(Fewbody can handle all these cases with ease since it 
is general in N). However, for simplicity we currently 
break triples into a binary and a single star. We do this 
by allowing the outer member of the triple to just barely 
escape to infinity, with the inner binary shrinking its or- 
bit to conserve energy in the process. For simplicity, both 
the single star and binary are given the systemic velocity 
of the original triple. 

Finally, we must discuss one more detail. Since the 
binary interactions performed with Fewbody are done 
in a vacuum — in other words, there are no stars other 
than the ones in the interaction to prevent members 
of the small-A'^ system from making arbitrarily large 
excursions — it happens that some binary interactions 
leave extremely wide binaries (sometimes as wide as the 
cluster itself) as their outcome products. Clearly this is 
an unphysical situation, as no binary should ever become 
larger than the inter-particle separation in a comparable- 
mass cluster. We therefore break these pathologically 
wide binaries at the end of each timestep. Our criterion 
for breaking the binaries is that they have an orbital ve- 
locity that is roughly smaller than the local velocity dis- 
persion, since it is this boundary in phase space and not 
the hard-soft boundary that determines binary lifetime 
in a cluster (Fregeau et al. 2006). We break the binaries 
if the orbital speed of their lightest member is less than 
Xhsa{R), where cr(i?) is the velocity dispersion at posi- 
tion R in the cluster, and Xy^s is a parameter which we 
take to be 0.7. We set < 1 as a safety measure, to 
ensure that no long-lived binaries arc erroneously broken. 

2.5. Stellar Evolution 

For the mass-radius relationship we adopt an approx- 
imate piece-wise fit to the more detailed one used by 
Freitag et al. (2006b) for Z = IQ-^: 



R{M) 



M/Mq < 0.1 
0.1 < M/Mq < 1 
1 < M/Mq < 120 ■ 
M/Mq > 120 

(12) 

The first two pieces (up to M = 1 Mq) are based on 
Chabricr & Baraffe (2000). The next piece is based on 
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Schallcr ct al. (1992), and the last is from Bond ct al. 
(1984). The fit is a rough approximation to the more 
exact relationships given in the references listed, but suf- 
fices for our purposes since we use it only for determining 
limits on the properties of the initial binary population 
(as described below). 

2.6. Timestep Evaluation 

The timestep in the code should be chosen small 
enough to resolve the relevant physics (two-body re- 
laxation, collisions, binary interactions, etc.), but not 
smaller than necessary. We take the timestep to be the 
minimum of the characteristic timescales for the different 
physical processes. 

For two-body relaxation, we use the standard expres- 
sion for the characteristic timescale for two species of 
particles undergoing relaxation to deflect each other by 
an angle ffmax (Freitag & Benz 2001): 



7r/2 32 \n{jN)G^n{Mi + A'hf 



(13) 



where Wrei is the relative speed of the two species, n is 
the local number density of stars, and Mi are the masses 
of each species of star. The standard expression for the 
relaxation time comes from setting ^max = 7r/2 (Binncy 
& Tremainc 1987). We evaluate eq. (13) by a local sliding 
average as 



7r/2 32 ln(7iV)G2n((A/i + M2Y) 



(14) 



yielding r,.ci as a function of radial position in the cluster. 
We take the minimum value of Tj-ei for the calculation of 
the timestep. The minimum most often occurs at the 
center of the cluster, where the density is the highest. 
However, it can sometimes happen for clusters with wide 
mass spectra that the minimum occurs away from the 
center, due to a massive star in a sea of lighter stars. 
We adopt ^max = 1 for all simulations presented in this 
paper, which we find to be a good compromise between 
accuracy and computational speed. The distribution of 
scattering angles in a typical timestep has a very long 
tail at large 9, so most super-encounters have a much 
smaller scattering angle than 0max- 

For strong interactions, we evaluate the timescale for 
each pair of particles neighboring in radius to undergo 
a strong interaction, by performing an "ncr?;" estimate. 
This can be written 

T,7/„„g = ~ / ^^Vid3v2/(V1)/(V2)|V2 - VilS-.trong , 

Q5) 

where v,; is the velocity of star i, and / is the velocity 
distribution function. Assuming a MaxwcUian velocity 
distribution for stars "1" and "2", the result is 



-1 

strong 



A^/irnria ( 1 



GM 
2rpcr2 



(16) 



where a is the one-dimensional velocity dispersion, and 
we have substituted eq. (9) for S'strong- For collisions, we 
plug in ScoW to find 

T-J, = 16V^n,i„g,e(i??)<T [l + , (17) 



where we explicitly show which quantities we average, 
and risingic is the number density of single stars. Similarly 
for binary-binary interactions: 

and binary-single interactions: 

T,-i = 4^/¥n.ingieX2,(a')^ [l + ^1^^) ' ^^^^ 

where nbin is the number density of binaries, and a is the 
binary semimajor axis. 

2.7. Initial Conditions 

For our initial cluster models we use both isolated 
Plummer models and tidally-truncated King models of 
varying concentration. Our prescription for tidal mass 
loss is described in detail in Joshi et al. (2001). For 
models with a mass spectrum or binaries we assume no 
primordial mass segregation of the heavier components. 
We use a binary fraction from to 1 , with the binary 
fraction defined as fb = Nb/{Ns + Ni,), where Ng is the 
number of single stars in the cluster, Nb is the number 
of binaries, and N = Ng + Nb. 

In assigning the initial properties of the binary popu- 
lation, we start with a cluster of only single stars. We 
create each binary by randomly choosing a cluster star 
to be the primary member of the binary, and assigning 
the secondary mass using a flat distribution for the bi- 
nary mass ratio q {dP/dq (x 1), truncated at the low end 
so that the mass of the secondary is not lower than the 
minimum of the initial mass function. With the masses 
of both binary members set, the remaining binary prop- 
erties are set according to one of two different schemes. 
The first is the scheme that has traditionally been used 
in numerical modeling of dense stellar systems, in which 
the binary binding energy Eb is distributed uniformly in 
the logarithm [dP/dEb (x E^^), with fixed upper and 
lower limits. As shown in Table 1, we take as limits on 
the binding energy a few kT on the low end to several 
hundred kT on the high end, where kT is the thermal en- 
ergy in the cluster core, evaluated as ^{mv^) « ("^)cr^, 
where is the one-dimensional velocity dispersion in 
the core. With the semimajor axis set by the binding 
energy, the eccentricity e is set according to the ther- 
mal distribution {dP/de = 2e). The second scheme is a 
slightly modified version of the first in which the limits 
on the binding energy and eccentricity are set in a more 
physical way. The binding energy is still distributed as 
dP/dEb (X E~^, but with the upper limit set to the bind- 
ing energy at a semimajor axis of 5(i?i + R2), where Ri 
arc the stellar radii. The lower limit on the binding en- 
ergy is set to that of a binary whose lightest member has 
orbital speed Xhs (vvci) , where Xhs is as defined in section 
2.4.2, and (urei) is the locally averaged relative velocity 
between objects, taken to be 4CT/37r, where a is the local 
three-dimensional velocity dispersion. The eccentricity 
is set according to the thermal distribution with an up- 
per limit set by the pericenter distance. Note that when 
physical limits on Eb are used, the resulting cluster sim- 
ulation is no longer scalable in its length, mass, and time 
units, since adopting stellar radii sets the physical scale 
of the system. 
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3. EXAMPLE RESULTS AND COMPARISONS 

Having verified that our code properly treats two-body 
relaxation and correctly samples the interaction rates 
for strong interactions, we now use it study the evolu- 
tion of more realistic clusters. We first consider clusters 
of equal-mass stars with primordial binary populations. 
Note again that for all simulations presented in this pa- 
per physical stellar collisions were turned off. The focus 
of this section is on presenting a few illustrative results in 
detail, and comparing the results of our newly modified 
code with those of other codes, as well as the previous 
version of our code. 

Fig. 9 shows the evolution of an isolated Plummer 
model with N = 10^ stars, and an initial 3% binary frac- 
tion (model pLnle5_fb0.03 in Table 1). The top panel 
shows Alb, the total mass in binaries bound to the cluster 
(solid line), and M, the total mass of the cluster (dashed 
line) as a function of time, relative to their initial val- 
ues. The middle panel shows i?bb, the cumulative energy 
generated in binary-binary interactions (solid line) and 
£^bs, the cumulative energy generated in binary-single 
interactions (dashed line) relative to |i^c,o|, the absolute 
value of the cluster's initial mechanical energy. The bot- 
tom panel shows the evolution of Vc, the cluster core 
radius (solid line), ni,b, the half-mass radius of the bina- 
ries (dashed line), and rh,s, the half- mass radius of single 
stars (dot-dashed line). The evolution of this model is 
typical of a cluster with primordial binaries. The core 
initially shrinks until the central density increases to the 
point at which energy generation in binary interactions 
is sufiicient to prevent the core from collapsing. The bi- 
naries steadily gain binding energy in the subsequent, 
long-lived binary burning phase. They thus suffer pro- 
gressively larger kinetic recoil from binary interactions, 
with the result that eventually the half-mass radius of bi- 
naries overtakes the single star half-mass radius. Eventu- 
ally the binary population is sufficiently depleted in the 
core that the core collapses, by which point a large frac- 
tion of the initial binary population has been depleted 
(^ 90%). The binaries are lost either by being disrupted 
or ejected in binary interactions. From the middle panel 
it is clear that a significant fraction of energy can be re- 
leased via binary interactions, in this case of order the 
initial cluster mechanical energy. 

Fig. 10 shows the evolution of the virial ratio and the 
total cluster energy (cluster mechanical energy plus bi- 
nary binding energy) for the same model. The virial 
ratio is conserved to within statistical fluctuations for 
the duration of the calculation, suggesting that the code 
is yielding accurate results. The total cluster energy is 
also conserved relatively well throughout the calculation, 
at the level of a part in ~ 10^, with the largest jump 
in the energy occurring during the deep core collapse 
phase. The high degree of conservation of energy and of 
the virial ratio for this model is typical of all models we 
present in this paper. 

3.1. Comparison with Theory 

For the first quantitative test of our results for the 
quasi steady-state binary burning phase, we appeal to the 
semi-analytical model of Vesperini & Chernoff (1994). 
Their model combines the results of binary scattering 
experiments with an analytical prescription for energy 




Fig. 9. — Evolution of an isolated Plummer model with N = 10'' 
stars, and an initial 3% binary fraction (model pl_nlc5_fb0.03 in 
Table 1). The top panel shows Mi,, the total mass in binaries 
bound to the cluster (solid line), and M, the total mass of the 
cluster (dashed line) as a function of time, relative to their ini- 
tial values. The middle panel shows -Ebbi the cumulative energy 
generated in binary-binary interactions (solid line) and -Ebs, the 
cumulative energy generated in binary-single interactions (dashed 
line) relative to |i?c,oli the absolute value of the cluster's initial 
mechanical energy. The bottom panel shows the evolution of rc, 
the cluster core radius (solid line), bi the half-mass radius of the 
binaries (dashed line), and rb,a, the half-mass radius of single stars 
(dot-dashed line). Time is in initial half-mass relaxation times. 




Fig. 10. — Evolution of the virial ratio (top panel), and total 
cluster energy (bottom panel) for model pLnle5_fb0.03. The en- 
ergy plotted here includes not only the mechanical energy of the 
cluster, but also the binding energy of the binaries. Note the range 
on the y axis for the energy plot. In this run energy was conserved 
to a part in ~ 10*. 

balance in the core, yielding a relationship among the 
core binary fraction, the properties of the binary popu- 
lation, and the ratio r^/rii: 

r, ^ 0.1872 nhsM^ - M + ^^bh^l (£\ 

ru logio(7^) {l + ^bf \vh) VlOy ' 

(20) 

where is the core binary fraction, ^bs and /ibb are 
coefficients representing the energy generation rates in 
binary-single and binary-binary interactions for a given 
set of binary properties, Vc and vu are the core and half- 
mass velocity dispersions, and F parameterizes the ex- 
pansion rate of the core in terms of the half-mass relax- 
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ation time. Hcggic ct al. (2006) have shown that eq. (20) 
agrees well with the results of A^-body simulations of 
N = 4096 clusters with primordial binaries, with the 
values r = 9.4 and Vc/vh = v^, which are close to the 
canonical values of F = 11.5 and Vc/vh = adopted 
in Vesperini & Chernoff (1994). However, they find that 
the dependence of r^/r/j on A^ is steeper than eq. (20), 
in the sense that clusters with A^ > 10'^ have a systemat- 
ically smaller value of rc/vh than eq. (20) predicts, with 
the opposite true for clusters with fewer stars. Extrap- 
olating from the numerical results in Fig. 18 of Heggie 
et al. (2006), it appears that for N ~ 10*^ eq. (20) over- 
estimates Tc/rh by a factor of « 2 for the values F = 10 
and Vc/vh = V^- 

To test the agreement between our code and eq. (20), 
we have performed several cluster simulations of varying 
initial binary fraction and measured the core radius and 
binary fraction after core stabilization. The details of the 
models (pl_n3e5_fb0.01_kt through pLn3e5_fb0.60_kt) are 
given in Table 1. Fig. 11 shows Vc/rh vs. core binary frac- 
tion (/);, for the simulations, compared with eq. (20). Each 
set of points (shown as circles and triangles, alternat- 
ing) represents a simulation with a different initial binary 
fraction. Each point is determined by averaging over a 
time window of width At = t^^ from the point of core 
stabilization until several tj-h later. The solid line shows 
the theoretical model with the standard values F = 10 
and Vc/vh = \/2 for Ei, = 10-100 fcT. The agreement 
between our code and the semi-analytical model is quite 
satisfactory for </);, < 0.5. Above this value the two differ 
by up to a factor of ~ 2, with our code yielding rc/rt 
values smaller than predicted by theory. We believe the 
apparent discrepancy is due to the destruction in binary- 
binary interactions of the wider binaries in the initial 
binary distribution, the rate of which increases with in- 
creasing binary fraction. For reference, the theoretical 
curve for a narrower range of binary binding energies 
{Eb = 25-100 fcT) is shown in the dashed line. 

Finally, we mention the issue of post-collapse, 
gravothermal core oscillations. Gravothcrmal oscillations 
are driven by the gravothermal instability (essentially the 
negative heat capacity of gravitational systems, which 
causes heat to flow from cold to hot in a runaway fash- 
ion). At any given time during deep core collapse, the 
temperature profile is a monotonically decreasing func- 
tion of r, with the central temperature steadily increasing 
as a function of time. At some point during the collapse, 
a binary scattering interaction occurs, producing a small 
amount of energy in the core and cooling it, creating a 
temperature inversion. Once the temperature inversion 
is established, the gravothermal instability takes over and 
drives the subsequent core expansion. There are several 
hallmarks of gravothermal oscillations. One is, of course, 
a temperature inversion in the cluster core at the point 
when the core begins to rebound. Another is that the 
expansion phase is not driven by energy generation in 
binaries. Yet another is the presence of loops in a phase 
space diagram of the cluster core properties (Makino 
1996; Heggie et al. 2006). In Paper HI we demonstrated 
the gravothermal nature of the core oscillations produced 
by our code by displaying the temperature inversion in 
the core. As in Paper HI, cluster models evolved with our 
newly modified code undergo core oscillations after core 




Fig. 11. — The quantity rc/r^ vs. core binary fraction (pi, for 
simulations pl_n3e5_fb0.01_kt through pl_n3e5_fb0.60_kt, compared 
with the semi-analytical model of Vesperini & Chernoff (1994). 
Each set of points (shown as circles and triangles, alternating) rep- 
resents a simulation with a different initial binary fraction (from 
left to right): = 0.01, 0.02, 0.04, 0.08, 0.15, 0.3, and 0.6. Each 
point is determined by averaging over a time window of width 
At = tjij from the point of core stabilization until several f^jj later. 
The solid line shows the theoretical model with the standard values 
r = 10 and Vc/v^ = \f2 for E^, = 10-100 fcT, while the dashed line 
shows the theoretical model with i?;, = 25-100 kT to reflect the 
depletion of the widest binaries due to binary-binary interactions. 

collapse. Preliminary tests show loops in phase space, 
but with the wrong directional sense, implying that the 
oscillations we now see are not gravothcrmal in nature. 
Since we are concerned only with pre-deep core collapse 
evolution in this paper, we postpone for future work a 
more detailed analysis of the core oscillations produced 
by our newly modified code. 

3.2. Comparison with Previous Numerical Work 

The amount of numerical work reported on clusters 
with primordial binaries has grown considerably since 
the work of Gao et al. (1991), who used a multimass 
Fokker-Planck code coupled with a recipes-based treat- 
ment of binary interactions. In Paper HI, we used a 
Monte Carlo code coupled with recipes for binary inter- 
actions. Giersz & Spurzcm (2003) used a hybrid code, 
which treated single stars via a gas dynamical method 
and binaries via Monte Carlo, and performed direct nu- 
merical integration of binary interactions. Heggie et al. 
(2006) and Trenti et al. (2006b) performed direct A^-body 
simulations. And as described in the preceding sections, 
in this paper we use a Monte Carlo method coupled with 
direct numerical integration of binary interactions. In 
this section wc compare the results from our code with 
those from all the methods just listed. 

A standard model has emerged in the business of simu- 
lating the evolution of clusters containing primordial bi- 
naries. The "Gao, ct al." model is an isolated Plummcr 
model with equal-mass stars and 10% binaries, with the 
binary binding energy distributed uniformly in the loga- 
rithm from 3 to 400A:T, and with the binary eccentricity 
distributed according to the thermal distribution. Since 
this model has been treated by all the previous work cited 
above, we use it as a basis for comparison. Fig. 12 shows 
the evolution of this model (T4 in Table 1). Fig. 13 shows 
the evolution of the ratio Tc/r/j. Although we did not in- 
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tegrate our model to deep core collapse, which occurs at 
t w 50frh in Gao et al. (1991) and t w ISOtrh in Heggie 
et al. (2006), there are still ample results for comparison. 
Comparing with Fig. 1 of Gao et al. (1991), the timescale 
of the initial core contraction we find is nearly identical 
at w llirh- This is the same timescale found by Heggie 
et al. (2006) in their Fig. 2 and Giersz & Spurzem (2003) 
in their Fig. 4, although in the latter case the comparison 
is less meaningful since their results show no long-lived 
binary burning phase. It is also similar to the result of 
Paper 111 in Fig. 4, although in that paper we found a 
slightly shorter timescale. 

Turning now to the size of the core radius during the 
binary burning phase, we find Tc ~ 0.07 in A^-body units. 
The model of Gao et al. (1991) shows a smaller core size, 
with Tc « 0.05 in iV-body units (note that the unit of 
length in their plot is 37r/16 in TV-body units). This is 
somewhat expected, since the code of Gao et al. (1991) 
predicts a smaller core than the semi-analytical theory 
(Heggie et al. 2006). The recipes-based model of Paper 
HI shows a much larger core, with Vc ~ 0.2. This is 
not surprising, since recipes tend to overestimate the en- 
ergy generation rate in binary interactions (Fregeau et al. 
2003; Giersz & Spurzem 2003; Fregeau et al. 2005). Al- 
though the model of Giersz & Spurzem (2003) does not 
show a long-lived binary burning phase, there is a subtle 
hint of a short-lived one with ~ 0.2. This is also much 
larger than the value we find. Comparing our Fig. 13 
with Fig. 2 of Heggie et al. (2006), we see qualitatively 
very similar behavior in the two calculations, with rc/vh 
abruptly slowing its contraction at t w ll^rh and steadily 
decreasing thereafter. At the start of the binary-burning 
phase, Heggie et al. (2006) find rc/r^ ~ 0.06, just as 
we do. If the scaling of rc/vh with N in Heggie et al. 
(2006) can be extrapolated out to iV ~ 10^ as described 
in section 3.1, then this agreement is somewhat surpris- 
ing. In terms of the structural parameters during the 
binary burning phase (which, as the longest-lived evo- 
lutionary phase in the life of a cluster is the most ob- 
servationally relevant), as well as the timescale to reach 
it, our code appears to agree well with A^-body, slightly 
less well with the Fokker-Planck code, and much less well 



with the two other approximate codes. 

We can also compare the mass in binaries (or, equiva- 
lently, number) as a function of time. At the somewhat 
arbitrary time of 40trh, we find M^/Mbfi k, 0.35. Gao 
et al. (1991) find a value of 0.41, while Heggie et al. 
(2006) find a value of 0.36. In Paper III we found a value 
of 0.35. All these methods appear to agree very well 
on the average rate that binaries are lost (either ejected 
or disrupted in binary interactions) up to 40irh7 although 
the agreement with Paper HI is probably fortuitous given 
the disagreement in the structural parameters. 

Comparing an isolated cluster model tests in combina- 
tion how well our code treats two-body relaxation, mass 
segregation, and the various aspects of binary burning. 
Having shown that our code agrees well with the "exact" 
method of iV-body, one can be fairly confident that our 
code treats these processes relatively accurately. How- 
ever, since all dense star clusters are tidally truncated 
to some degree, it is still useful to consider how well our 
code treats the physics of tidal stripping. As described 
in Joshi et al. (2001), we use a cutoff radius criterion to 
determine whether a star has been stripped from the clus- 
ter. This is in contrast to the more accurate technique of 
including the tidal field in the equations of motion (which 
is not possible for the Monte Carlo method). Trenti et al. 
(2006b) have compared the two tidal stripping methods 
by performing A^-body simulations of iV ~ 10'* clusters 
with primordial binaries, and have found that models us- 
ing a tidal cutoff tend to survive longer before disrupting 
(up to a factor of ^ 2, although it's not clear how this fac- 
tor scales with N) , and have a larger core radius (again, 
up to a factor of ~ 2), than models with a tidal field. 
We will keep this in mind when presenting results below. 
For now we compare two current simulations with the 
tidal cutoff models of Trenti et al. (2006b), and with the 
previous version of our code. Figs. 14 and 15 show the 
evolution of models w3_nle5_fb0.1 and w7_nle5_fb0.1, re- 
spectively. Comparing the first with Fig. 13 in Paper III, 
we see that although both models reach disruption, our 
new results are qualitatively different in behavior, yield- 
ing what appears to be a binary burning phase from 9irh 
to disruption, while the old result shows no such phase. 
Aside from this difference, the quantitative evolution of 
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Fig. 14.— Evolution of model THH3. Quantities plotted are the FiG. 15.— Evolution of model THH7. Quantities plotted are the 

same as in Fig. 9. same as in Fig. 9. 



the structural radii, the total cluster and binary mass, 
and the disruption timcscalc appear to be very similar 
between the two models. Comparing with Fig. 16 of 
Trenti et al. (2006b), we find similar qualitative behav- 
ior, with a hint of a binary burning phase in their results 
starting at roughly the same time, similar evolution of 
the structural radii, and both models resulting in dis- 
ruption at ~ 12-14trh. Our model seems to predict a 
smaller core radius, by a factor of ~ 2. Note that Trenti 
ct al. (2006b) use the Spitzer definition of the core radius, 
which they find with their code to yield a value ^ 20% 
larger than the density-weighted averaging method. 

Moving on to the Wq = 7 model, we find similar qual- 
itative behavior but shorter disruption times than both 
Fig. 10 in Paper III and Trenti et al. (2006b) in their 
Fig. 17. The recipes-based model of Paper III predicts 
a much larger core radius in the binary burning phase 
than our current model. This is due to the fact that the 
recipes used in that calculation tend to overestimate the 
rate of energy generation in binary burning. The agree- 
ment with Trenti et al. (2006b) is better, with our core 
radius being only ~ 20% smaller than the iV-body re- 
sult. Looking at the evolution of M and Aff,, our model 
appears to exhibit very similar behavior to the iV-body 
model. 

In general, our new code yields much improved agree- 
ment with the A'^-body results reported in the literature, 
for both isolated and tidally-truncated models, but yields 
some results that differ significantly from other approxi- 
mate methods. Notably, the direct integration of binary 
interactions reduces their energy generation rate relative 
to the simple recipes used in Paper III, and yields smaller 
core radii. This was evident to some degree in Paper III, 
in which the discrepancy was illustrated for binary-single 
interactions. Below we compare our new results with ob- 
servations and discuss the implications. 

4. RESULTS AND COMPARISON WITH OBSERVATIONS 

In addition to displaying the initial conditions for all 
models simulated for this paper. Table 1 gives several 
important measured quantities for each simulation. The 
first is the core stabilization time, tcs, at which the core 
radius stabilizes (i.e., the start of the binary-burning 
phase) after the initial contraction or expansion. Note 




Fig. 16. — Evolution of model wll_nle5_fb0.3. Quantities plot- 
ted are the same as in Fig. 9. 
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Fig. 17. — Evolution of rc/rh (dashed line), total cluster binary 
fraction (dot-dashed line) , core binary fraction /]-, ^ (dotted line) , 
and concentration c (solid line) for model wll_nle5_fb0.3. 
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Fig. 18. — Evolution of c and rc/r/^ for Wo = 3, 7, and 11 models 
with /i, = 0.1 (solid lines) and = 0.3 (dashed lines). Note the 
"universal" behavior, with the initially more centrally concentrated 
models (Wq = 11) expanding and the less centrally concentrated 
models (Wo = 3 and 7) contracting to common values of rc/ri^ and 
c. 




5 10 15 20 25 30 35 40 

Fig. 19. — Evolution of the binary binding energy distribution 
for model T5. The distribution starts off uniform in logSj, from 
3 to 400 kT. As the cluster evolves the central velocity disper- 
sion increases, moving the hard-soft boundary to larger Ei, and 
destroying the wider binaries. Near core collapse (~ 40trh), only a 
small collection of relatively tight binaries {Ei, ~ 300 kT) remains. 



that some authors denote this as the core collapse time 
when discussing clusters with primordial binaries. The 
second is the time of the first deep core collapse, tec, 
which, for models with binaries, represents the time at 
which the binary population is nearly depleted in the 
core. The next is the disruption time, idis, for models 
which are tidally truncated. The next is ^c/r/j, the ra- 
tio of the core to half-mass radius averaged over 1 tj-h 
after t^s- Finally, c ~ log]^Q(7'f /r^) is the concentration 
parameter averaged over the same time period. The evo- 
lution of several of the models in the tabic has been 



shown graphically in the preceding figures. The behav- 
ior of the remaining models is similar to the ones already 
shown, with the relevant timescales and structural pa- 
rameters appropriately modified. The only exceptions 
arc the high concentration King models, which undergo 
an initial phase of core expansion (instead of contraction) 
due to their initially very dense state. Fig. 16 displays 
the evolution of such a model, model wll_nle5_fb0.3. For 
reference we have plotted the time evolution of the mea- 
sured quantities Vc/rh and c in Fig. 17, along with the 
the total cluster binary fraction /(, and the core binary 
fraction /b,c- 

Fig. 18 shows the evolution of c and rc/vh for Wq = 3, 
7, and 11 models with /& = 0.1 (sohd lines) and fb = 0.3 
(dashed lines). Of note is the "universal" behavior dis- 
played by the models, in which the initially more cen- 
trally concentrated models (Wq = H) expand and the 
less centrally concentrated models (Wq = 3 and 7) con- 
tract to common values of Vc/rh and c. This is the same 
behavior found previously in the literature, in Paper III 
and Trenti ct al. (2006b) in their Fig. 8, although in 
the latter case they find systematically larger values of 
rc/rh and smaller values of c (as described above). For 
reference we also show in Fig. 19 the evolution of the 
binary binding energy distribution for model T5. The 
softest binaries begin to be destroyed (by ejection or dis- 
ruption in binary interactions) on the approach to core 
stabilization. As the cluster evolves the central velocity 
dispersion increases, moving the hard-soft boundary to 
larger E}, and destroying the wider binaries. Near core 
collapse (^ 40<rh), only a small collection of relatively 
tight binaries {E^ ~ 300 kT) remains. The evolution of 
the binary population in E^-r space is qualitatively sim- 
ilar to what's shown in Fig. 15 of Paper III, with the 
softer binaries first being destroyed in the core, and the 
harder binaries becoming even harder and more centrally 
concentrated with time. 

There arc several trends apparent in the data presented 
in Table 1. For fixed initial cluster structure, as fb is 
increased from zero, tcs decreases from the value of tec 
at fb = 0, reaches a minimum around /{, ~ 0.2, then 
increases back to approximately its fb = value for 
/(, = !. This dip is due to the fact that in the initial 
core contraction or expansion phase, binaries act pri- 
marily as a second, heavier star species, hastening clus- 
ter energy transport via mass segregation. Looking at 
the deep core collapse times, another clear trend is that 
tec increases dramatically as fb is increased. This trend 
is most striking in the isolated Plummer models, with 
merely fb ~ 0.03 being enough to double the deep core 
collapse time, and fb = ^ increasing it by a factor of 
at least 7 (model T4). The physical explanation is. of 
course, that the lifetime of the binary-burning phase in- 
creases with the amount of fuel available, in analogy to 
hydrogen-burning in main-sequence stars. Looking more 
carefully at the tidally-truncated King models, we see 
that the presence of binaries tends to drive these clus- 
ters to complete tidal disruption. The minimum of <dis 
occurs somewhere in the range 0.1 < fb < 1. This is 
not surprising, since the maximum of rc/rh occurs at 
fb ~ 0.4 (see Fig. 11), implying that the cluster is the 
most distended for this value of fb- Moving now to the 
cluster structural parameters, we see that there is rela- 
tively little variation in rc/rh and c over the wide range 
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of cluster initial profiles and binary fractions considered, 
with Tc/rh ranging from ~ 0.05 for Plummer or Wq ~ 7 
King models with low fb, to ~ 0.1 for larger Wq King 
models with ff, ^ 1. The concentration parameter shows 
a similarly small amount of variation, peaking at ~ 2.1 
for Wo ~ 7 and low binary fraction, and falling to ^ 1.7 
for larger Wq and larger binary fraction. As expected, a 
larger /f, leads to a smaller c, since the core radius gen- 
erally increases with binary fraction. Another trend is 
evident when comparing single-mass models with models 
incorporating a more realistic mass spectrum (Salpeter 
from 0.2-1.2 Mq): models with a mass spectrum tend to 
have a slightly smaller Vc/rh, and show more variation of 
c with Wq. 

Although our cluster evolution models are rather sim- 
plified (since they do not include single- or binary-star 
stellar evolution, or collisions), it is still useful to com- 
pare the predicted structural parameters with observa- 
tions. In Paper III we compared rc/rh and c from sim- 
ulations using the previous version of our Monte Carlo 
code (which includes recipes for binary interactions in- 
stead of dynamical integrations), with observations for 
the Galactic globular clusters. There we found promis- 
ing agreement, with rc/rji from the simulations falling 
generally in the low rc/rh region of the observed dis- 
tribution for non-core collapsed clusters, which extends 
from ~ 0.1 to ~ 1 with a peak at ^ 0.5. And similarly 
for c, falling in the high c region of the observed dis- 
tribution for non-core collapsed clusters, which extends 
from ^ 0.5 to ~ 2.5 with a peak at ~ 1.5. However, 
as described above, we find systematically smaller values 
for Tc/rh than in Paper III by a factor of up to ^ 10, 
as well as systematically larger values for c by ~ 0.5. 
The improved treatment of binary interactions in our 
code has shifted our predictions for rc/rh and c outside 
the observed ranges for non-core collapsed clusters, now 
yielding agreement only with the roughly 10% of Galac- 
tic globular clusters that are classified observationally as 
core collapsed. 

A globular cluster of comparable mass stars with /f, > 
0.03 viewed at a random time during its life has a greater 
than 50% chance of being found in the quasi-steady state 
binary burning phase. When one considers that globu- 
lars are likely born with significant binary fractions (Hut 
et al. 1992; Ivanova et al. 2005), and that we are currently 
observing the Galactic globular clusters at a late stage 
in their evolution, the vast majority of observed clusters 
should currently be in the binary-burning phase. The 
most obvious interpretation of the observational data is 
that most Galactic globular clusters are currently in the 
binary-burning phase, and the roughly 10% classified as 
core-collapsed are within a small time window around a 
deep core collapse phase. The disagreement between sim- 
ulations and observations for the structural parameters 
of the non core-collapsed clusters, then, suggests one of 
at least two possibilities: 1) the Galactic globular clus- 
ters do not start within the volume in parameter space of 
initial conditions we have considered here; or 2) there are 
additional physical processes at work in clusters, yield- 
ing larger cores. If clusters are born with Wq ^11 and 
/b < 1, extrapolation of our results suggests that simu- 
lations may then agree with the observations of non-core 
collapsed clusters, implying possibility (1) may be cor- 
rect. However, a cluster with Wo > 11 is likely to have 



a short enough central relaxation time that a runaway 
stellar collision will occur, creating an intermediate mass 
black hole (IMBH) early in the cluster's lifetime (Fre- 
itag et al. 2006a). This is intriguing, since clusters with 
central IMBHs typically have significantly larger values 
of Tc/r/i than clusters without (Trcnti et al. 2006a). We 
have not included in our simulations any form of stellar 
evolution (for single stars or binaries), or physical stellar 
collisions. Single star evolution tends to heat a cluster 
early in its lifetime via wind-driven mass loss and su- 
pernovae explosions, causing the cluster and its core to 
expand. However, this effect is most pronounced only 
early in the lifetime of a star cluster (< IGyr). The 
effects of binary stellar evolution are less obvious, since 
it is a rather complicated process. However, the simu- 
lations of Ivanova et al. (2005) suggest that the binary 
fraction in the core quickly drops to relatively small val- 
ues (< 20%), and that binary stellar evolution tends to 
destroy tight binaries. The net result is likely to be a 
smaller equilibrium value of rc/rh (Vespcrini & Chcrnoff 
1994) than with no binary stellar evolution, suggesting 
that possibility (2) is less likely. A refined quantitative 
study combining the results of Ivanova et al. (2005) and 
Vespcrini & Chernoff (1994) would more clearly elucidate 
this effect. 

The effect of direct single-single star collisions in young 
dense clusters is to dissipate orbital energy and drive core 
collapse (Freitag et al. 2006a). For clusters in which stel- 
lar merger products have had time to evolve and lose 
mass through accelerated stellar evolution, the net re- 
sult may be to heat the core (Lee 1987). The degree 
to which this process operates in Galactic globular clus- 
ters is unclear, however. The effect of stellar collisions 
during binary interactions is generally to reduce the effi- 
ciency of binary burning (Hut & Inagaki 1985; McMillan 
1986; Goodman & Hcrnquist 1991). This is because a 
merger product resulting from a star-star collision typ- 
ically has significantly more internal energy (potential 
and rotational) than the sum of the merging stars' inter- 
nal energies. Thus when a collision occurs during a bi- 
nary dynamical interaction, effectively some of the bind- 
ing energy of the binary is converted into stellar binding 
energy, decreasing the efficiency of binary burning. In 
other words, energy that had previously been available 
to be converted into kinetic energy through binary inter- 
actions is no longer available, tied up in the stars. The 
result should be a decreased value of rc/rh in the binary- 
burning phase relative to the case of point-mass binary 
dynamics, implying that possibility (2) is less likely. As 
with binary stellar evolution, however, a more detailed 
simulation which studies the effects of collisions in an 
evolving model should be performed to quantify the ef- 
fect. 

5. SUMMARY AND CONCLUSIONS 

In this paper we have described our new Monte Carlo 
evolution code and used it to perform a large set of clus- 
ter evolution simulations, which wc compared with pre- 
vious results in the literature, as well as observations of 
Galactic globular clusters. 

In section 2 we described our new code in detail, in- 
cluding the implementation of direct integration of bi- 
nary scattering interactions and star-star physical colli- 
sions, and the fundamental modifications we have made 
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to the core Monte Carlo method. We performed sev- 
eral test calculations with the code, finding that it re- 
produces well several standard results. It yields a core 
collapse time of w 18 ^rh for an isolated Plummer model, 
in good agreement with results in the literature. It also 
produces an r"^ "^ density profile during the late stages 
of core collapse, in good agreement with the theoretical 
expectation. We also performed comparisons of clusters 
with increasingly wide mass spectra with A'^-body, find- 
ing that for moderately wide mass spectra (1 to 10 Mq) 
the agreement with A^-body is satisfactory, but for very 
wide mass spectra (0.2 to 120 Mq) the agreement is not 
as good. In particular, for such wide mass spectra, our 
Monte Carlo code tends to overestimate the mass segre- 
gation timescale at early times, and underestimate it at 
later times. The sense of the disagreement is the same as 
found by Freitag et al. (2006b) with their Monte Carlo 
code. 

In section 3 we displayed a few example results and 
compared with theory and previous numerical calcula- 
tions in the literature. We found that the code con- 
serves energy well over the long timescales of our runs. 
We compared our predicted values of rc/rh during the 
quasi-steady state binary burning phase with the semi- 
analytical work of Vesperini & ChernofF (1994), finding 
good agreement. We compared our results with previous 
numerical results in the literature for isolated and tidally- 
truncated cluster models, finding excellent agreement 
with iV-body calculations. There are much larger dis- 
crepancies with the other approximate methods (Fokker- 
Planck and other Monte Carlo codes), which is to be ex- 
pected, since most used recipes for binary interactions, 
which are known to overestimate the energy generation 
rate in binaries. 

In section 4 we surveyed the results from all our sim- 
ulations, and compared with observations. Our mod- 
els cover a large range in parameter space, using Plum- 
mer and King models with Wq = 3 to 11 for the initial 
cluster profile, and with initial binary fractions from 
to 1. The resulting structural parameters in the binary 



burning phase span a remarkably small range, with rc/rh 
varying from 0.03 to 0.12, and c varying from 1.7 to 2.4. 
Our results for these structural parameters are distinctly 
different from the results found with the previous ver- 
sion of our code (which used recipes for binary interac- 
tions), with rc/rh now smaller than what we found in 
Paper III by a factor of up to ^ 10, and with c larger 
by ^ 0.5. Although our new results agree much bet- 
ter with A'-body calculations, they unfortunately agree 
much less well than in Paper III with the observations. 
The disagreement implies one of at least two possibilities. 
It may be that the initial conditions for Galactic globu- 
lar clusters are outside the range of initial conditions we 
have sampled in this work. Extrapolation of our results 
suggests that clusters with H^o H and /& ^ 1 may 
match the observations. This is intriguing since a cluster 
with Wq > 11 will likely form an IMBH early in its life- 
time via a collisional runaway, and clusters with central 
IMBHs are known to have larger values of rc/rh (Freitag 
et al. 2006a; Trenti et al. 2006a). Alternatively, stellar 
evolution and collisions, which are not included in the 
simulations in this paper, could possibly explain the dis- 
agreement. However, it appears that the effect of these 
processes should act in the opposite sense of ameliorating 
the disagreement with observations. More detailed sim- 
ulations including the effects of single- and binary-star 
evolution and physical collisions should be performed to 
test this. 
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TABLE 1 

Parameters and measured quantities for all model simulations presented in this paper. 



name profile tnb/pc /(Af/Afo) ft f(Ei,/kT) ics/trh tcc/irh tdis/trh rc/rh c 



Tl 


5 X 10'"^ 


Plum. 


1.02 


(X S{M - 1) 











17.6 








T2 


10'5 


Plum. 


0.58 


Kroup., 0.1-10 











0.54 








T3 


1.25 X 10'' 


Plum. 


0.60 


Salp., 0.2-120 











0.067 








T4 


10^ 


Plum. 


7.49 


(X 5{M - 1) 


0.1 




3-400 


11 


> 128 




0.06 




T5 


10^ 


Plum. 


7.33 


(X 5{M - 1) 


0.03 




3-400 


14 


42 




0.06 




THH3 


10^ 


VKo =3 


7.49 


(X 5{M - 1) 


0.1 




3-400 


9.4 


11.8 


> 11.8 


0.07 


1.8 


THH7 


10^ 


Wo =7 


7.49 


(X 5{M - 1) 


0.1 


(xi?-i 


3-400 


6.7 


29.8 


> 29.8 


0.08 


1.9 


pLn3e5.fb0.01.kt 


3 X 10^ 


Plum. 


5.44 


(X 5{M - 1) 


0.01 


(X E-\ 


10-100 


16.3 


> 23 




0.04 




pLn3e5.fb0.02_kt 


3 X 10^ 


Plum. 


5.46 


(X 5{M - 1) 


0.02 


(X E-^, 


10-100 


14.2 


> 23 




0.06 




pLn3e5_fb0.04_kt 


3 X 10-"^ 


Plum. 


5.49 


(X 5{M - 1) 


0.04 


° 1 

(X E-^, 
1 


10-100 


12.5 


> 23 
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Note. — Here A*^ is the number of total cluster objects (single stars and binaries), the profile is either a Plummer model or a King model 
with the specified Wo, r'NB is the unit of length in the simulation, f{M) is the initial mass function, is the initial binary fraction, f{E)j) is 
the distribution of binary binding energy, tcs is the time at which the core radius stabilizes (i.e., the start of the binary- burning phase) after 
the initial contraction or expansion (note that some authors denote this as the core collapse time when discussing clusters with primordial 
binaries), tec is the time of the first deep core collapse, t^i^ is the time at which the cluster disrupts due to tidal stripping, rc/ri^ is the ratio 
of the core to half-mass radius averaged over Itjij after tcs, and c = logj^Q (rt/rc) is the concentration parameter averaged over the same time 
period. Quantities are omitted when the physical state they describe is never reached (or can never be reached) during the simulation. Note 
that those models with = or those with fcT-based limits on Ei, have one degree of freedom in their scaling, while those with physical 
limits on the binary population cannot be rescaled. 



